# Replicate paper empirics (descriptive statistics and analysis)

# Setup ----
# setwd("") # copy-paste here your working directory if you do not work in the R project

# library(maps)
library(showtext)
library(kableExtra)
library(modelsummary)
library(interflex)
library(tidyverse)
library(tidylog, warn.conflicts = FALSE)

## Plot defaults ----
# setup font of plots to be Times New Roman using showtext
font_add("Times New Roman", "times new roman.ttf") # use this if you're on a UNIX system (Mac/Linux)
# font_add("Times New Roman", "times.ttf") # use this if you're on a Windows system
showtext_auto()

# ggplot theme default:
theme_set(theme(legend.position = "bottom", legend.direction = 'horizontal',
                panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                panel.background = element_rect(fill = "transparent", color = NA),
                plot.background = element_rect(fill = "transparent", color = NA),
                legend.background = element_rect(fill = "transparent", color = NA),
                legend.box.background = element_rect(fill = "transparent", color = NA),
                legend.key = element_rect(fill = "transparent", color = NA), 
                axis.line = element_line(colour = "black", linewidth = 0.3),
                text = element_text(size = 16, family = "Times New Roman"),
                legend.text = element_text(size = 12),
                legend.title = element_text(size = 12, hjust = 0.5)))

# Table defaults ----
options("modelsummary_format_numeric_latex" = "plain")

cf_map <- c("post_treat:treatment" = "AR6 $\\times$ WGI-III",
            "no_delegates_tot" = "Delegation size",
            "post_treat" = "AR6",
            "treatment" = "WGI-III",
            "(Intercept)" = "(Intercept)")

# Import data ----
partic <- read_rds("data/data_participation.rds") 
engage <- read_rds("data/data_engagement.rds")

# Descriptives ----
## P. 2 statements: country participation in AR5 vs AR6 ----
partic %>%
  filter(no_delegates_tot > 0) %>%
  select(AR, country) %>%
  distinct() %>%
  group_by(AR) %>%
  summarise(n = n())
# number of countries participating per AR

partic %>%
  group_by(country, AR) %>%
  summarise(no_delegates_tot = sum(no_delegates_tot)) %>%
  filter(country %in% country[AR == "AR5" & no_delegates_tot > 0] &
           country %in% country[AR == "AR6" & no_delegates_tot == 0]) %>%
  select(country) %>%
  distinct() %>%
  nrow()
# 18 countries with at least one delegate in a session of AR5 did not send any delegate to any session of AR6

partic %>%
  group_by(country, AR) %>%
  summarise(no_delegates_tot = sum(no_delegates_tot)) %>%
  filter(country %in% country[AR == "AR5" & no_delegates_tot == 0] &
           country %in% country[AR == "AR6" & no_delegates_tot > 0]) %>%
  select(country) %>%
  distinct() %>%
  nrow()
# 31 countries with no delegate in any session of AR5 did send at least one delegate to a session of AR6

## Supplementary Material A list of countries participating to both AR5 and AR6 ----
d <- partic %>%
  filter(no_delegates_tot > 0) %>%
  group_by(country, AR) %>%
  summarise(no_delegates_tot = sum(no_delegates_tot)) %>%
  ungroup() %>%
  pivot_wider(id_cols = "country", names_from = "AR", values_from = "no_delegates_tot") %>%
  mutate(country_AR5 = case_when(!is.na(AR5) ~ country),
         country_AR6 = case_when(!is.na(AR6) ~ country)) %>%
  select(-AR6, -AR5) %>%
  mutate(order = !is.na(country_AR5) & !is.na(country_AR6)) %>%
  arrange(-order, country_AR5, country_AR6) %>%
  select(-country, -order) %>%
  rename("AR5 Participant Countries" = "country_AR5",
         "AR6 Participant Countries" = "country_AR6")

d %>%
  filter(!is.na(`AR5 Participant Countries`) & 
           !is.na(`AR6 Participant Countries`)) %>%
  select(1) %>%
  pull() %>%
  str_replace_all("\\&", "\\\\&") %>%
  paste0(collapse = ", ") %>%
  write("tables/si_list.tex")

## Table SI1 ----
t <- d %>%
  filter(!is.na(`AR5 Participant Countries`) & 
           is.na(`AR6 Participant Countries`)) %>%
  select(1) %>%
  rename("Countries participating only to AR5" = 1) %>%
  mutate(n = row_number()) %>%
  full_join(d %>%
              filter(is.na(`AR5 Participant Countries`) & 
                       !is.na(`AR6 Participant Countries`)) %>%
              select(2) %>%
              rename("Countries participating only to AR6" = 1) %>%
              mutate(n = row_number())) %>%
  select(-n) %>%
  mutate(N = row_number() %>% sprintf("%1.0f", .),
         `Countries participating only to AR5` = case_when(is.na(`Countries participating only to AR5`) ~ "",
                                                           TRUE ~ `Countries participating only to AR5`)) %>%
  relocate(N) %>%
  datasummary_df(title = "List of countries that sent at least one delegate to any IPCC negotiation only in AR5 or only in AR6 \\label{tab:countries}",
                 escape = TRUE,
                 output = "latex") %>%
  kable_styling(latex_options = c("hold_position"),
                font_size = 10, full_width = FALSE)

# export:
t %>% save_kable("tables/tab_si1.tex")

## Figure 1 ----
partic %>%
  filter(working_group != "SYN") %>%
  group_by(working_group, AR, treatment) %>%
  summarise(no_delegates_tot = mean(no_delegates_tot, na.rm = TRUE)) %>% # these are the individual WGs' averages
  # add the pooled WGs averages:
  bind_rows(partic %>%
              filter(working_group != "SYN") %>%
              group_by(treatment, AR) %>%
              summarise(no_delegates_tot = mean(no_delegates_tot, na.rm = TRUE)) %>%
              mutate(working_group = "Pooled Working Groups")) %>%
  # add the control (SYN) averages:
  bind_rows(partic %>%
              filter(working_group == "SYN") %>%
              group_by(treatment, AR) %>%
              summarise(no_delegates_tot = mean(no_delegates_tot, na.rm = TRUE)) %>%
              mutate(working_group = "I")) %>%
  bind_rows(partic %>%
              filter(working_group == "SYN") %>%
              group_by(treatment, AR) %>%
              summarise(no_delegates_tot = mean(no_delegates_tot, na.rm = TRUE)) %>%
              mutate(working_group = "II")
  ) %>%
  bind_rows(partic %>%
              filter(working_group == "SYN") %>%
              group_by(treatment, AR) %>%
              summarise(no_delegates_tot = mean(no_delegates_tot, na.rm = TRUE)) %>%
              mutate(working_group = "III")
  ) %>%
  bind_rows(partic %>%
              filter(working_group == "SYN") %>%
              group_by(treatment, AR) %>%
              summarise(no_delegates_tot = mean(no_delegates_tot, na.rm = TRUE)) %>%
              mutate(working_group = "Pooled Working Groups")
  ) %>%
  mutate(AR = str_extract(AR, '\\d') %>% as.numeric(),
         working_group = factor(working_group, levels = c("Pooled Working Groups", "I", "II", "III"),
                                labels = c("Pooled Working Groups", "Working Group I", "Working Group II", "Working Group III"))) %>%
  ggplot(aes(x = AR, y = no_delegates_tot, linetype = as.factor(treatment), shape = as.factor(treatment))) +
  geom_point() +
  geom_line() +
  lemon::facet_rep_wrap("working_group", repeat.tick.labels = "all", ncol=4) +
  xlab("Assessment Report") + ylab("Delegation size (average)") +
  coord_cartesian(ylim = c(0, 4)) +
  scale_x_continuous(breaks = c(5, 6), labels = function(x) paste0("AR", x)) +
  scale_linetype_manual("", breaks = c(1, 0), labels = c("Working Groups\n(online in AR6)", "Synthesis Report\n(in-person always)"),
                        values = c("solid", "dotted")) +
  scale_shape_manual("", breaks = c(1, 0), labels = c("Working Groups\n(online in AR6)", "Synthesis Report\n(in-person always)"),
                     values = c(19, 1)) +
  guides(linetype = guide_legend(keywidth = unit(1, "cm")))
ggsave("plots/fig_1.pdf", height = 5, width = 9)

## Pp. 7-8 statements on participation ----
tab <- partic %>%
  filter(working_group != "SYN") %>%
  group_by(country, AR) %>%
  summarise(avg_no_delegates_tot = mean(no_delegates_tot)) %>%
  pivot_wider(id_cols = "country", names_from = "AR", values_from = "avg_no_delegates_tot")

tab <- tab %>%
  filter(!(AR5 == 0 & AR6 == 0)) %>%
  ungroup()

nrow(tab) # 156 countries that sent at least one delegate to either AR5 or AR6 working group sessions

tab %>%
  filter(AR6 > AR5) %>%
  nrow() # 93 increased their average delegation size in AR6

round(93/156*100) # 60%

tab %>%
  filter(AR6 < AR5) %>%
  nrow() # 52 increased their average delegation size in AR6

round(52/156*100) # 33%

tab %>%
  filter(AR6 == AR5) %>%
  nrow() # 11 increased their average delegation size in AR6

round(11/156*100) # 7%

tab %>%
  filter(AR6 > AR5) %>%
  mutate(diff = AR6 - AR5) %>%
  select(diff) %>%
  pull() %>%
  mean() %>%
  round(digit = 1) # for countries that increased delegation size, average increase is +3.6 

tab %>%
  filter(AR6 < AR5) %>%
  mutate(diff = AR6 - AR5) %>%
  select(diff) %>%
  pull() %>%
  mean() %>%
  round(digit = 1) # for countries that decreased delegation size, average decrease is -0.8 

tab %>%
  filter(AR6 < AR5) %>%
  mutate(diff = AR6 - AR5) %>%
  filter(diff < -1) 
# only 6 countries reduced their participation by more than 1 delegate. Japan
# is the one with largest reduction (hosting WGII in Yokohama in AR5)

tab %>%
  filter(AR5 != 0 & AR6 == 0) %>%
  nrow() # 31 countries in WG sessions of AR5 but not in those of AR6

tab %>%
  filter(AR6 != 0 & AR5 == 0) %>%
  nrow() # 25 countries in WG sessions of AR6 but not in those of AR5

tab %>%
  filter(country %in% c("Antigua & Barbuda", "Samoa", "St. Kitts & Nevis", "Vanuatu"))

tab %>%
  mutate(diff = AR6 - AR5,
         perc = diff/AR5) %>%
  arrange(-AR6) %>%
  print(n = Inf)

tab %>%
  mutate(diff = AR6 - AR5,
         rate = diff/AR5) %>%
  filter(AR5 != 0) %>%
  arrange(-rate) %>%
  print(n = Inf)

## Figure SI1 ----
partic %>%
  group_by(country_code, working_group, country) %>%
  summarise(no_delegates_AR5 = no_delegates_tot[AR == "AR5"],
            no_delegates_AR6 = no_delegates_tot[AR == "AR6"]) %>%
  ungroup() %>%
  mutate(diff = no_delegates_AR6 - no_delegates_AR5,
         sign = case_when(diff > 0 ~ "positive",
                          diff == 0 ~ "same",
                          diff < 0 ~ "negative"),
         working_group = factor(working_group, levels = c("SYN", "I", "II", "III"),
                                labels = c("Synthesis Report", "Working Group I", "Working Group II", "Working Group III")),
         country = tidytext::reorder_within(country, -diff, working_group),
         pos_label = case_when(sign == "positive" ~ no_delegates_AR6 + 5,
                               sign == "same" ~ no_delegates_AR6 + 3,
                               sign == "negative" ~ no_delegates_AR5 + 3)) %>%
  ggplot(aes(x = country)) +
  geom_segment(aes(xend = country, y = no_delegates_AR5, yend = no_delegates_AR6, col = sign),
               arrow = arrow(length = unit(0.05, "inches"), type = "closed")) +
  geom_text(aes(y = pos_label, label = country_code, col = sign), size = 2, show.legend = FALSE,
            family = "Times New Roman") +
  lemon::facet_rep_wrap("working_group", scales = "free_x", repeat.tick.labels = "all") +
  # coord_flip() +
  tidytext::scale_x_reordered() +
  expand_limits(x = c(-5, 166)) +
  xlab("Country ISO-2 code") + ylab("Number of delegates in AR5 and AR6") +
  scale_color_discrete("", breaks = c("positive", "same", "negative"), labels = c("Increased size", "Same size", "Decreased size")) +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank())
ggsave("plots/fig_si1.pdf", width = 12, height = 6)

## Table 1 ----
did1 <- lm(no_delegates_tot ~ post_treat * treatment, 
           data = partic)

did2 <- lm(no_delegates_tot ~ post_treat * treatment, 
           data = partic %>%
             filter(working_group %in% c("I", "SYN")))

did3 <- lm(no_delegates_tot ~ post_treat * treatment, 
           data = partic %>%
             filter(working_group %in% c("II", "SYN")))

did4 <- lm(no_delegates_tot ~ post_treat * treatment, 
           data = partic %>%
             filter(working_group %in% c("III", "SYN")))

tab_did1 <- modelsummary(models = list("Pooled" = did1, "WGI" = did2,  "WGII" = did3, "WGIII" = did4),
                         stars = TRUE,
                         coef_map = cf_map,
                         gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",
                         output = "latex",
                         title = "Effect of virtual negotiations on average delegation size (DID estimator) \\label{tab:models_DID}",
                         escape = FALSE) %>%
  kable_styling(latex_options = c("hold_position"),
                font_size=10, full_width = FALSE) %>% 
  add_header_above(c(" ", "Model 1" = 1, "Model 2" = 1, "Model 3" = 1, "Model 4" = 1)) %>%
  footnote(general_title="Notes:",
           footnote_as_chunk = FALSE,
           threeparttable = TRUE,
           escape = FALSE,
           general = c("Outcome: Delegation size. Standard errors in parentheses. Model~1 pools data across WGs; models 2$-$4 show estimates for WGs separately."))

save_kable(tab_did1, file = "tables/tab_1.tex")

## Figure 2 ----
facet_dat <- engage %>%
  filter(working_group != "SYN") %>%
  group_by(working_group) %>% 
  slice_max(order_by = no_delegates_tot, n = 30) %>%
  left_join(partic %>%
              filter(AR == "AR5" & working_group != "SYN") %>%
              select(country_code, working_group, no_delegates_tot) %>%
              rename("no_delegates_tot_AR5" = "no_delegates_tot"),
            by = c("country_code", "working_group")) %>%
  mutate(working_group = factor(working_group, levels = c("I", "II", "III"),
                                labels = c("Working Group I", "Working Group II", "Working Group III"))) %>%
  pivot_longer(cols = c("no_delegates_male", "no_delegates_female", "no_delegates_NA", "tot_interventions"),
               names_prefix = "no\\_delegates\\_",
               names_to = "gender",
               values_to = "count") %>%
  mutate(facet = case_when(gender == "tot_interventions" ~ "interventions",
                           TRUE ~ "delegation"),
         count = case_when(gender == "tot_interventions" ~ -count,
                           TRUE ~ count),
         gender = factor(gender, levels = c("female", "tot_interventions", "male", "NA")), 
         facet = factor(facet, levels = c("interventions", "delegation"),
                        labels = c("Engagement", "Size of delegation")))

p1 <- facet_dat %>%
  filter(working_group == "Working Group I") %>%
  bind_rows(data.frame(country = "...", no_delegates_tot = 0, gender = NA,
                       working_group = "Working Group I",
                       facet = c("Engagement", "Size of delegation"))) %>%
  ggplot(aes(x = reorder(country, no_delegates_tot))) +
  geom_histogram(aes(y = count, fill = gender), stat = "identity") +
  geom_point(data = . %>%
               filter(facet == "Size of delegation"),
             aes(y = no_delegates_tot_AR5, shape = "AR5 Size of delegation"), size = 2) +
  ggh4x::facet_nested(. ~ working_group + facet, scales = "free") +
  tidytext::scale_x_reordered() +
  coord_flip() +
  xlab("") + ylab("") +
  geom_hline(yintercept = 0) +
  scale_y_continuous(labels = function(x) ifelse(x < 0, abs(x), x)) +
  scale_fill_manual("", values = c("female" = scales::hue_pal()(3)[[1]], 
                                   "male" = scales::hue_pal()(3)[[3]],
                                   "NA" = "grey40", 
                                   "tot_interventions" = "grey"),
                    breaks = c("female", "male", "NA")) +
  scale_shape_manual("", values = 4) +
  ggh4x::facetted_pos_scales(y = list(facet == "Engagement" ~ scale_y_continuous(limits = c(-250, 0), 
                                                                                 labels = function(x) ifelse(x < 0, abs(x), x)),
                                      facet == "Size of delegation" ~ scale_y_continuous(limits = c(0, 30)))) +
  theme(panel.spacing = unit(-0.75, "lines"))

p2 <- facet_dat %>%
  filter(working_group == "Working Group II") %>%
  bind_rows(data.frame(country = "...", no_delegates_tot = 0, gender = NA,
                       working_group = "Working Group II",
                       facet = c("Engagement", "Size of delegation"))) %>%
  ggplot(aes(x = reorder(country, no_delegates_tot))) +
  geom_histogram(aes(y = count, fill = gender), stat = "identity") +
  geom_point(data = . %>%
               filter(facet == "Size of delegation"),
             aes(y = no_delegates_tot_AR5, shape = "AR5 Size of delegation"), size = 2) +
  ggh4x::facet_nested(. ~ working_group + facet, scales = "free") +
  tidytext::scale_x_reordered() +
  coord_flip() +
  xlab("") + ylab("") +
  geom_hline(yintercept = 0) +
  scale_y_continuous(labels = function(x) ifelse(x < 0, abs(x), x)) +
  scale_fill_manual("", values = c("female" = scales::hue_pal()(3)[[1]], 
                                   "male" = scales::hue_pal()(3)[[3]],
                                   "NA" = "grey40", 
                                   "tot_interventions" = "grey"),
                    breaks = c("female", "male", "NA")) +
  scale_shape_manual("", values = 4) +
  ggh4x::facetted_pos_scales(y = list(facet == "Engagement" ~ scale_y_continuous(limits = c(-250, 0), 
                                                                                 labels = function(x) ifelse(x < 0, abs(x), x)),
                                      facet == "Size of delegation" ~ scale_y_continuous(limits = c(0, 30)))) +
  theme(panel.spacing = unit(-0.74, "lines"))

p3 <- facet_dat %>%
  filter(working_group == "Working Group III") %>%
  bind_rows(data.frame(country = "...", no_delegates_tot = 0, gender = NA,
                       working_group = "Working Group III",
                       facet = c("Engagement", "Size of delegation"))) %>%
  ggplot(aes(x = reorder(country, no_delegates_tot))) +
  geom_histogram(aes(y = count, fill = gender), stat = "identity") +
  geom_point(data = . %>%
               filter(facet == "Size of delegation"),
             aes(y = no_delegates_tot_AR5, shape = "AR5 Size of delegation"), size = 2) +
  ggh4x::facet_nested(. ~ working_group + facet, scales = "free") +
  tidytext::scale_x_reordered() +
  coord_flip() +
  xlab("") + ylab("") +
  geom_hline(yintercept = 0) +
  scale_y_continuous(labels = function(x) ifelse(x < 0, abs(x), x)) +
  scale_fill_manual("", values = c("female" = scales::hue_pal()(3)[[1]], 
                                   "male" = scales::hue_pal()(3)[[3]],
                                   "NA" = "grey40", 
                                   "tot_interventions" = "grey"),
                    breaks = c("female", "male", "NA")) +
  scale_shape_manual("", values = 4) +
  ggh4x::facetted_pos_scales(y = list(facet == "Engagement" ~ scale_y_continuous(limits = c(-250, 0), 
                                                                                 labels = function(x) ifelse(x < 0, abs(x), x)),
                                      facet == "Size of delegation" ~ scale_y_continuous(limits = c(0, 30)))) +
  theme(panel.spacing = unit(-0.74, "lines"))

ggpubr::ggarrange(p1, p2, p3, nrow = 1, common.legend = TRUE, legend = "bottom",
                  align = "h", widths = c(0.36, 0.32, 0.32))
dev.copy("plots/fig_2.pdf", height = 8, width = 15, device = pdf)
dev.off()

## Pp. 10-11 statements on participation and engagement ----
partic %>%
  filter(country == "Japan")

partic %>%
  filter(country == "Saudi Arabia")

partic %>%
  filter(country == "Germany")

partic %>%
  filter(country == "China")

df <- engage %>%
  filter(working_group != "SYN") %>%
  select(country, working_group, no_delegates_tot) %>%
  group_by(working_group) %>% 
  slice_max(order_by = no_delegates_tot, n = 30)

t.test(df$no_delegates_tot[df$working_group == "I"], df$no_delegates_tot[df$working_group == "II"]) # 9.6 vs 12.8, p = 0.025
t.test(df$no_delegates_tot[df$working_group == "I"], df$no_delegates_tot[df$working_group == "III"]) # 9.6 vs 13.8, p = 0.006
t.test(df$no_delegates_tot[df$working_group == "II"], df$no_delegates_tot[df$working_group == "III"]) # 12.8 vs 13.8, p = 0.534

engage %>%
  filter(working_group != "SYN") %>%
  group_by(working_group) %>%
  mutate(tot = sum(tot_interventions)) %>%
  slice_max(order_by = no_delegates_tot, n = 30) %>%
  select(country, working_group, tot_interventions, tot) %>%
  group_by(working_group) %>%
  summarise(tot_top30 = sum(tot_interventions),
            tot = unique(tot)) %>%
  mutate(perc = 100*tot_top30 / tot)

engage %>%
  filter(working_group != "SYN") %>%
  group_by(working_group) %>%
  mutate(tot = sum(tot_interventions),
         IN_SA_US = case_when(country %in% c("India", "Saudi Arabia", "United States") ~ 1,
                              TRUE ~ 0)) %>%
  select(country, working_group, tot_interventions, tot, IN_SA_US) %>%
  group_by(working_group, IN_SA_US) %>%
  summarise(tot_group = sum(tot_interventions),
            tot = unique(tot)) %>%
  mutate(perc = 100 * tot_group / tot)

engage %>%
  group_by(country, working_group) %>%
  summarise(tot = sum(tot_interventions)) %>%
  arrange(working_group, -tot) %>%
  print(n = Inf)

(sum(engage$tot_interventions[engage$working_group == "III"]) /
    sum(engage$tot_interventions[engage$working_group == "I"])) %>%
  round()
# 2 times

(sum(engage$tot_interventions[engage$working_group == "III"]) /
    sum(engage$tot_interventions[engage$working_group == "II"])) %>%
  round()
# 2 times

## Figure 3 ----
engage %>%
  filter(working_group == "II" | working_group == "III") %>%
  # comment out the following line to keep the full list of delegations:
  group_by(working_group) %>% slice_max(order_by = share_out_of_office, n = 30) %>%
  bind_rows(data.frame(country = "...", share_out_of_office = 0,
                       working_group = c("II", "III"))) %>%
  mutate(working_group = factor(working_group, levels = c("II", "III"),
                                labels = c("Working Group II", "Working Group III")),
         country = tidytext::reorder_within(country, share_out_of_office, working_group)) %>%
  ggplot(aes(x = country, y = share_out_of_office)) +
  geom_histogram(stat = "identity") +
  facet_wrap("working_group", scales = "free_y") +
  tidytext::scale_x_reordered() +
  coord_flip() +
  xlab("") + ylab("Share of negotiation hours outside of 9am-5pm office hours")
ggsave("plots/fig_3.pdf", device = "pdf", height = 6, width = 10)

engage %>%
  filter(working_group %in% c("II", "III")) %>%
  filter(country %in% c("Vanuatu", "Australia", "South Korea", "Japan", "Cook Islands",
                        "Samoa", "New Zealand", "Kiribati")) %>%
  select(country, working_group, tot_out_of_office, tot_in_office) %>%
  mutate(perc = 100*tot_out_of_office / (tot_out_of_office + tot_in_office)) %>%
  arrange(-perc)

## Figure 4 ----
int1 <- interflex(estimator = "binning",
                  data = engage  %>%
                    mutate(l_tot_interventions = log(tot_interventions+1)) %>% as.data.frame(),
                  Y = "l_tot_interventions", D = "no_delegates_tot", X = "share_out_of_office",
                  treat.type = "continuous", na.rm = TRUE, figure = FALSE)
rescale <- 0.3
shift1 <- 0.3
p1 <- int1[["est.lin"]][[1]] %>%
  as.data.frame() %>%
  as_tibble() %>%
  mutate(across(.cols = c("ME", "lower CI(95%)", "upper CI(95%)"),
                .fns = ~.x+shift1),
         working_group = "Pooled Working Groups") %>%
  ggplot() +
  geom_line(aes(x = X, y = ME)) +
  geom_ribbon(aes(x = X, ymin = `lower CI(95%)`, ymax = `upper CI(95%)`),
              alpha = I(.3)) +
  geom_pointrange(data = int1[["est.bin"]][[1]] %>% 
                    as.data.frame() %>%
                    as_tibble() %>%
                    mutate(across(.cols = c("coef", "CI.lower", "CI.upper"),
                                  .fns = ~.x+shift1),
                           working_group = "Pooled Working Groups"),
                  mapping = aes(x = x0, y = coef, ymin = CI.lower, ymax = CI.upper)) +
  geom_hline(yintercept = 0+shift1, linetype = "dashed") +
  geom_histogram(data = engage %>%
                   mutate(working_group = "Pooled Working Groups"),
                 aes(x = share_out_of_office,
                     y = rescale*((after_stat(count))/sum(after_stat(count))))) +
  scale_y_continuous(labels = function(x) x - shift1) +
  coord_cartesian(ylim = c(0, 0.8)) +
  facet_wrap("working_group") +
  xlab("Share of negotiation hours out of office") +
  ylab("Marginal effect of delegation size on engagement")

int2 <- interflex(estimator = "binning",
                  data = engage %>% 
                    filter(working_group == "II") %>%
                    mutate(l_tot_interventions = log(tot_interventions+1)) %>% as.data.frame(),
                  Y = "l_tot_interventions", D = "no_delegates_tot", X = "share_out_of_office",
                  treat.type = "continuous", na.rm = TRUE, figure = FALSE)

shift2 <- 0.3
p2 <- int2[["est.lin"]][[1]] %>%
  as.data.frame() %>%
  as_tibble() %>%
  mutate(across(.cols = c("ME", "lower CI(95%)", "upper CI(95%)"),
                .fns = ~.x+shift2),
         working_group = "Working Group II") %>%
  ggplot() +
  geom_line(aes(x = X, y = ME)) +
  geom_ribbon(aes(x = X, ymin = `lower CI(95%)`, ymax = `upper CI(95%)`),
              alpha = I(.3)) +
  geom_pointrange(data = int2[["est.bin"]][[1]] %>% 
                    as.data.frame() %>%
                    as_tibble() %>%
                    mutate(across(.cols = c("coef", "CI.lower", "CI.upper"),
                                  .fns = ~.x+shift2),
                           working_group = "Working Group II"),
                  mapping = aes(x = x0, y = coef, ymin = CI.lower, ymax = CI.upper)) +
  geom_hline(yintercept = 0+shift2, linetype = "dashed") +
  geom_histogram(data = engage %>%
                   filter(working_group == "II") %>%
                   mutate(working_group = "Working Group II"),
                 aes(x = share_out_of_office,
                     y = rescale*((after_stat(count))/sum(after_stat(count))))) +
  scale_y_continuous(labels = function(x) x - shift2) +
  coord_cartesian(ylim = c(0, 0.8)) +
  facet_wrap("working_group") +
  xlab("Share of negotiation hours outside of 9am-5pm office hours") +
  ylab("Marginal effect of delegation size on engagement")

int3 <- interflex(estimator = "binning",
                  data = engage %>% 
                    filter(working_group == "III") %>%
                    mutate(l_tot_interventions = log(tot_interventions+1)) %>% as.data.frame(),
                  Y = "l_tot_interventions", D = "no_delegates_tot", X = "share_out_of_office",
                  treat.type = "continuous", na.rm = TRUE, figure = FALSE)

shift3 <- 0.3
p3 <- int3[["est.lin"]][[1]] %>%
  as.data.frame() %>%
  as_tibble() %>%
  mutate(across(.cols = c("ME", "lower CI(95%)", "upper CI(95%)"),
                .fns = ~.x+shift3),
         working_group = "Working Group III") %>%
  ggplot() +
  geom_line(aes(x = X, y = ME)) +
  geom_ribbon(aes(x = X, ymin = `lower CI(95%)`, ymax = `upper CI(95%)`),
              alpha = I(.3)) +
  geom_pointrange(data = int3[["est.bin"]][[1]] %>% 
                    as.data.frame() %>%
                    as_tibble() %>%
                    mutate(across(.cols = c("coef", "CI.lower", "CI.upper"),
                                  .fns = ~.x+shift3),
                           working_group = "Working Group III"),
                  mapping = aes(x = x0, y = coef, ymin = CI.lower, ymax = CI.upper)) +
  geom_hline(yintercept = 0+shift3, linetype = "dashed") +
  geom_histogram(data = engage %>%
                   filter(working_group == "III") %>%
                   mutate(working_group = "Working Group III"),
                 aes(x = share_out_of_office,
                     y = rescale*((after_stat(count))/sum(after_stat(count))))) +
  scale_y_continuous(labels = function(x) x - shift3) +
  coord_cartesian(ylim = c(0, 0.8)) +
  facet_wrap("working_group") +
  xlab("Share of negotiation hours out of office") +
  ylab("Marginal effect of delegation size on engagement")

ggpubr::ggarrange(p1 + xlab(""), 
                  p2 + ylab(""), 
                  p3 + xlab("") + ylab(""), nrow = 1)
ggsave("plots/fig_4.pdf", width = 10, height = 5)

## Export session info ----
sessioninfo::session_info(to_file = "session.log")

#====# The End #====#